function show_vfield(Vx, Vy, intv, roi)
%SHOW_VFIELD Shows a velocity field 
%
%   SHOW_VFIELD(Vx, Vy, intv);
%
%       Shows a velocity field. Here, Vx and Vy are the maps of the
%       x-component values and the y-component values. 
%
%       intv is the interval between adjacent grid points where the
%       arrows are shown, which can be either a scalar or a pair in 
%       form of [intv_x, intv_y].
%
%   SHOW_VFIELD(Vx, Vy, intv, roi);
%
%       Only shows within the region of interests. (roi is a logical mask).
%

% Created by Dahua Lin, on April 5, 2012

%% verify inputs 

if ~(isfloat(Vx) && isreal(Vx) && ismatrix(Vx))
    error('show_vfield:invalidarg', 'Vx should be a real matrix.');
end
if ~(isfloat(Vy) && isreal(Vy) && ismatrix(Vy))
    error('show_vfield:invalidarg', 'Vx should be a real matrix.');
end

[h, w] = size(Vx);
if ~isequal(size(Vy), [h, w])
    error('show_vfield:invalidarg', 'Vx and Vy should be of the same size.');
end

if isscalar(intv)
    intv_x = intv;
    intv_y = intv;
elseif numel(intv) == 2
    intv_x = intv(1);
    intv_y = intv(2);
else
    error('show_vfield:invalidarg', 'The form of intv is invalid.');
end

if nargin < 4
    use_roi = false;
else
    if ~(islogical(roi) && isequal(size(roi), [h w]))
        error('show_vfield:invalidarg', 'roi should be a h x w logical matrix.');
    end
    use_roi = true;
end


%% main

xv = 1 + round(rem(w-1, intv_x) / 2 : intv_x : w-1);
yv = 1 + round(rem(h-1, intv_y) / 2 : intv_y : h-1);

[xx, yy] = meshgrid(xv, yv);

idx = sub2ind([h, w], yy, xx);

if use_roi
    si = find(roi(idx));
    
    xx = xx(si);
    yy = yy(si);
    idx = idx(si);
end

vx = Vx(idx);
vy = Vy(idx);

quiver(xx, yy, vx, vy);





